Attribute VB_Name = "Options"
Option Base 1
Option Explicit
Public Function EuropeanOption(CallPut As String, S As Double, X As Double, r As Double, time As Double, sigma As Double) As Double

    If Left(CallPut, 1) = "c" Then
        EuropeanOption = EuropeanCall(S, X, r, time, sigma)
    ElseIf Left(CallPut, 1) = "p" Then
        EuropeanOption = EuropeanPut(S, X, r, time, sigma)
    Else
        'MsgBox "Invalid Option Type " & CallPut & " in EuropeanOption", vbOKOnly, "European Option"
    End If
    
End Function
Public Function EuropeanCall(S As Double, X As Double, r As Double, time As Double, sigma As Double) As Double
'   S - spot price
'   X - Strike (exercise) price,
'   r - risk-free interest rate
'   time - Fractions of year
'   sigma - Annualized percent volatility
    EuropeanCall = EuropeanCallP(S, X, r, time, sigma)
    
End Function
Private Function EuropeanCallP(S As Double, X As Double, r As Double, time As Double, sigma As Double) As Double

    Dim time_sqrt As Double, D1 As Double, D2 As Double
    time_sqrt = Sqr(time)
    D1 = (Log(S / X) + (r + 0.5 * sigma ^ 2) * time) / (sigma * time_sqrt)
    D2 = D1 - (sigma * time_sqrt)
    EuropeanCallP = S * Application.WorksheetFunction.NormSDist(D1) - X * Exp(-r * time) * Application.WorksheetFunction.NormSDist(D2)

End Function
Public Function EuropeanPut(S As Double, X As Double, r As Double, time As Double, sigma As Double) As Double
'   S - spot price
'   X - Strike (exercise) price,
'   r - risk-free interest rate
'   time - Fractions of year
'   sigma - Annualized percent volatility

    Dim time_sqrt As Double, D1 As Double, D2 As Double
    time_sqrt = Sqr(time)
    D1 = (Log(S / X) + (r + 0.5 * sigma ^ 2) * time) / (sigma * time_sqrt)
    D2 = D1 - (sigma * time_sqrt)
    EuropeanPut = X * Exp(-r * time) * Application.WorksheetFunction.NormSDist(-D2) - S * Application.WorksheetFunction.NormSDist(-D1)

End Function
Public Function ImpliedVolatility(S As Double, X As Double, r As Double, time As Double, option_price As Double) As Double

'   Implied Volatility Calculation using Newton's method:
'       Xi+1 = Xi - f(X)/f'(X)

'   S - spot price
'   X - Strike (exercise) price,
'   r - risk-free interest rate
'   time - Fractions of year
'   option_price - Option Value
  
    '   check for arbitrage violations:
    '   if price at almost zero volatility greater than price, return 0

    Const sigma_low = 0.00001
    Const MAX_ITERATIONS = 100
    Const ACCURACY = 0.0001
    Dim price As Double, time_sqrt As Double, sigma As Double, i As Integer, diff As Double, D1 As Double, vega As Double
    price = EuropeanCallP(S, X, r, time, sigma_low)
    If (price > option_price) Then
        ImpliedVolatility = 0#
        Exit Function
    End If
    time_sqrt = Sqr(time)
    sigma = (option_price / S) / (0.398 * time_sqrt) ' find initial value
    For i = 0 To MAX_ITERATIONS
        price = EuropeanCallP(S, X, r, time, sigma)
        diff = option_price - price
        If (Abs(diff) < ACCURACY) Then
            ImpliedVolatility = sigma
            Exit Function
        End If
        D1 = (Log(S / X) + r * time) / (sigma * time_sqrt) + 0.5 * sigma * time_sqrt
        vega = S * time_sqrt * (1# / 2.506628275) * Exp(-0.5 * D1)      ' This is the first derivative of NORMSDIST w.r.t. sigma
        sigma = sigma + diff / vega
    Next i
    ImpliedVolatility = -99#

End Function
Public Function KirksApproximation(CallPut As String, S1 As Double, S2 As Double, X As Double, r As Double, time As Double, sigma1 As Double, sigma2 As Double, rho As Double)
    
    If Left(CallPut, 1) = "c" Then
        KirksApproximation = KirksApproximationCall(S1, S2, X, r, time, sigma1, sigma2, rho)
    ElseIf Left(CallPut, 1) = "p" Then
        KirksApproximation = KirksApproximationPut(S1, S2, X, r, time, sigma1, sigma2, rho)
    Else
        'MsgBox "Invalid Option Type " & CallPut & " in EuropeanOption", vbOKOnly, "European Option"
    End If
    
End Function
Public Function KirksApproximationCall(S1 As Double, S2 As Double, X As Double, r As Double, time As Double, sigma1 As Double, sigma2 As Double, rho As Double)

'   S1 - spot price on Asset 1
'   S2 - spot price on Asset 2
'   X - Strike (exercise) price,
'   r - risk-free interest rate
'   time - Fractions of year
'   sigma1 - Annualized percent volatility of Asset 1
'   sigma2 - Annualized percent volatility of Asset 2
'   rho - Correlation in price returns of Asset 1 and Asset 2

    Dim SA As Double, sigma As Double, time_sqrt As Double, D1 As Double, D2 As Double
    SA = S1 / (S2 + X)
    sigma = Sqr(sigma1 ^ 2 + (sigma2 * S2 / (S2 + X)) ^ 2 - 2 * rho * sigma1 * sigma2 * S2 / (S2 + X))
    time_sqrt = Sqr(time)
    D1 = (Log(SA) + (0.5 * sigma ^ 2) * time) / (sigma * time_sqrt)
    D2 = D1 - (sigma * time_sqrt)
    KirksApproximationCall = (S2 + X) * Exp(-r * time) * (SA * Application.WorksheetFunction.NormSDist(D1) - Application.WorksheetFunction.NormSDist(D2))

End Function
Public Function KirksApproximationPut(S1 As Double, S2 As Double, X As Double, r As Double, time As Double, sigma1 As Double, sigma2 As Double, rho As Double)

'   S1 - spot price on Asset 1
'   S2 - spot price on Asset 2
'   X - Strike (exercise) price,
'   r - risk-free interest rate
'   time - Fractions of year
'   sigma1 - Annualized percent volatility of Asset 1
'   sigma2 - Annualized percent volatility of Asset 2
'   rho - Correlation in price returns of Asset 1 and Asset 2

    Dim SA As Double, sigma As Double, time_sqrt As Double, D1 As Double, D2 As Double
    SA = S1 / (S2 + X)
    sigma = Sqr(sigma1 ^ 2 + (sigma2 * S2 / (S2 + X)) ^ 2 - 2 * rho * sigma1 * sigma2 * S2 / (S2 + X))
    time_sqrt = Sqr(time)
    D1 = (Log(SA) + (0.5 * sigma ^ 2) * time) / (sigma * time_sqrt)
    D2 = D1 - (sigma * time_sqrt)
    KirksApproximationPut = (S2 + X) * Exp(-r * time) * (Application.WorksheetFunction.NormSDist(-D2) - SA * Application.WorksheetFunction.NormSDist(-D1))

End Function
Public Function SpreadOption(CallPut As String, S1 As Double, S2 As Double, X As Double, r As Double, time As Double, sigma1 As Double, sigma2 As Double, rho As Double)

'   Monte-Carlo Valuation of a spread option:
'           Call--Max(S1-S2-X,0)
'           Put--Max(S2-S1-X,0)

'   S1 - spot price on Asset 1
'   S2 - spot price on Asset 2
'   X - Strike (exercise) price,
'   r - risk-free interest rate
'   time - Fractions of year
'   sigma1 - Annualized percent volatility of Asset 1
'   sigma2 - Annualized percent volatility of Asset 2
'   rho - Correlation in price returns of Asset 1 and Asset 2

    Dim Z1 As Double, Z2 As Double, W1 As Double, W2 As Double
    Dim Drift1 As Double, Drift2 As Double, Vol1 As Double, Vol2 As Double
    Dim St1 As Double, St2 As Double
    Dim Sum As Double, antirho As Double
    Dim i As Integer
    
    Const SimCount As Integer = 10000
    Rnd -1
    Randomize 5566
    Sum = 0#
    antirho = Sqr(1 - rho)
    Drift1 = (r - 0.5 * sigma1 ^ 2) * time
    Drift2 = (r - 0.5 * sigma2 ^ 2) * time
    Vol1 = sigma1 * Sqr(time)
    Vol2 = sigma2 * Sqr(time)
    
    For i = 1 To SimCount
        '   Use antithetical sampling to preserve zero mean
        If i Mod 2 = 1 Then
            W1 = Z(Rnd())
            W2 = rho * W1 + antirho * Z(Rnd())
        Else
            W1 = -W1
            W2 = -W2
        End If
        St1 = S1 * Exp(Drift1 + Vol1 * W1)
        St2 = S2 * Exp(Drift2 + Vol1 * W2)
        If CallPut = "c" Then
            If (St1 > St2 + X) Then Sum = Sum + St1 - St2 - X
        Else
            If (St2 > St1 + X) Then Sum = Sum + St2 - St1 - X
        End If
        'Debug.Print St1 & "," & St2 & "," & W1 & "," & W2
    Next i
    SpreadOption = Sum / SimCount
    
End Function
Private Function Z(ByVal p As Double) As Double

    ' ALGORITHM AS 111, APPL.STATIST., VOL.26, 118-121, 1977.
    ' PRODUCES NORMAL DEVIATE CORRESPONDING TO LOWER TAIL AREA = P.
    Dim q As Double, r As Double
    Const Split As Double = 0.42, A0 As Double = 2.50662823884, A1 As Double = -18.61500062529, A2 As Double = 41.39119773534, A3 As Double = -25.44106049637
    Const B1 As Double = -8.4735109309, B2 As Double = 23.08336743743, B3 As Double = -21.06224101826, B4 As Double = 3.13082909833
    Const C0 As Double = -2.78718931138, C1 As Double = -2.29796479134, C2 As Double = 4.85014127135, C3 As Double = 2.32121276858
    Const D1 As Double = 3.54388924762, D2 As Double = 1.63706781897
    Const Zero = 0, One = 1#, Half = 0.5
    
    q = p - Half
    If Abs(q) <= Split Then
        '0.08 < P < 0.92
        r = q * q
        Z = q * (((A3 * r + A2) * r + A1) * r + A0) / ((((B4 * r + B3) * r + B2) * r + B1) * r + One)
        Exit Function
    End If
    
    'P < 0.08 OR P > 0.92, Set R = MIN(P,1-P)
    r = p
    If q > Zero Then r = One - p
    If r <= Zero Then
        Z = 0
        Exit Function
    End If
    r = Sqr(-Log(r))
    Z = (((C3 * r + C2) * r + C1) * r + C0) / ((D2 * r + D1) * r + One)
    If q < Zero Then Z = -Z

End Function
Private Function NormCDF(ByVal Z As Double, Optional Upper As Boolean = True) As Double
    '  Algorithm AS66 Applied Statistics (1973) vol.22, no.3
    
    '  Evaluates the tail area of the standardised normal curve
    '  from X to infinity if Upper is true or missing, or
    '  from minus infinity to X if Upper is false
    
    ' Local constants
    Const Zero = 0#, One = 1#, Half = 0.5, Con = 1.28
    Dim Y#, Test#
    Dim up As Boolean
    
    '*** machine dependent constants
    Const ltone = 7#, utzero = 18.66
    
    Const p = 0.398942280444, q = 0.39990348504, _
          r = 0.398942280385, A1 = 5.75885480458, _
          A2 = 2.62433121679, A3 = 5.92885724438, _
          B1 = -29.8213557807, B2 = 48.6959930692, _
          C1 = -0.000000038052, C2 = 0.000398064794, _
          C3 = -0.151679116635, C4 = 4.8385912808, _
          C5 = 0.742380924027, C6 = 3.99019417011, _
          D1 = 1.00000615302, D2 = 1.98615381364, _
          D3 = 5.29330324926, D4 = -15.1508972451, _
          D5 = 30.789933034
    
    up = Upper
    If Z < Zero Then
        up = Not up
        Z = -Z
    End If
    If Not (Z <= ltone Or up) Or Z > utzero Then
        NormCDF = Zero
        If up Then NormCDF = One - NormCDF
        Exit Function
    End If
    Y = Half * Z * Z
    If Z <= Con Then
        Test = (p - q * Y / (Y + A1 + B1 / (Y + A2 + B2 / (Y + A3))))
        NormCDF = Half - Z * Test
        If up Then NormCDF = One - NormCDF
        Exit Function
    End If
    
    'NormCDF = r * Exp(-y) / (z + c1 + d1 / (z + c2 + d2 / (z + c3 + d3 / (z + c4 + d4 / (z + c5 + d5 / (z + c6))))))
    Test = Z + C3 + D3 / (Z + C4 + D4 / (Z + C5 + D5 / (Z + C6)))
    NormCDF = r * Exp(-Y) / (Z + C1 + D1 / (Z + C2 + D2 / Test))
    If up Then NormCDF = One - NormCDF

End Function


